#include "quadgaussseidelsolver.h"
#include <QApplication>
#include <QDir>
#include <QProcess>
#include <cstdio>
#include <iostream>
#include <cmath>

using namespace std;
QuadGaussSeidelSolver::QuadGaussSeidelSolver(int _n,double mu1, double e1, QObject * parent):
    BaseGaussSeidelSolver(_n, mu1, e1, parent)/*, k((n+2)*(n+2))*/
{
    k = m*m;
    QObject *par;
    par = new QObject();
    QProcess proc(par);

    proc.setWorkingDirectory(qApp->applicationDirPath() + QDir::separator());
    proc.start("python coeffs.py");
    proc.waitForStarted();
    proc.write(QString("%1\n").arg(2).toStdString().c_str());
    proc.write(QString("%1\n").arg(n).toStdString().c_str());
    proc.write(QString("%1\n").arg(mu1).toStdString().c_str());
    proc.write(QString("%1\n").arg(e1).toStdString().c_str());
    //proc.closeWriteChannel();

    proc.waitForReadyRead();
    zone_count = QString(proc.readLine()).toInt();
    //cout << zone_count << endl;


    //coeffs_list = new Zone[zone_count];
    //coeffs_map = new int[k];

    for (int zone = 0; zone < zone_count; zone++)
    {
        Zone temp(2);
        for (int i = 0; i < 2; i++)
        {
            temp[i].resize(2);
            for (int j = 0; j < 2; j++)
            {
                temp[i][j].resize(7);
                for (int p=0; p < 7; p++)
                {
                    proc.waitForReadyRead();
                    temp[i][j][p] = QString(proc.readLine()).toDouble();
                    //coeffs_list[zone].c[i][j][p] = QString(proc.readLine()).toDouble();
                    //cout << coeffs_list[zone].c[i][j][p] << " ";
                    //cout << temp[i][j][p] << " ";
                }
                //cout << endl;
            }

        }
        coeffs_list.push_back(temp);
        //cout << endl;
    }
    // reading coefficients map
    for (int i = 0; i < k; i++)
    {
        proc.waitForReadyRead();
        coeffs_map.push_back(QString(proc.readLine()).toInt());
        //coeffs_map[i] =QString(proc.readLine()).toInt();
        //cout << coeffs_map[i] << endl;
    }

//    x = new double[2*k];


    x.resize(2*k);
    for  (int i = 0; i < k; i++)
    {
        if (coeffs_map[i] != -1)
        {
            //x[i] = ((double)random() / RAND_MAX - 0.5);
            //x[i+k] = ((double)random() / RAND_MAX - 0.5);
            x[i] = x[i+k] = 0.0;
        }
        else
        {
            x[i] = x[i+k] = 0.0;
        }
    }
    cout << QString(proc.readLine()).toStdString() << endl;
    proc.waitForFinished();
}

double QuadGaussSeidelSolver::iteration()
{
    double maxdx = 0;
    for (int i = 0; i < k; i++)
    {
        Zone  b;
        if (coeffs_map[i] == -1)
            continue;
        b = coeffs_list[coeffs_map[i]];
//        for (int ii = 0; ii < 2; ii++)
//        {
//            for (int j = 0; j < 2; j++)
//            {
//                for (int p = 0; p < 7; p++)
//                    b[ii][j][p] = coeffs_list[coeffs_map[i]][ii][j][p];
//            }
//        }

        double x_new = 0.0;
        double y_new = 0.0;
        double newdx = 0;
        int column = i % m;
        int row = i / m;
        //cout << "row = " << row << "; column  = " << column;
        if (row == 1)
        {
            x[i] = 0;
            x[i+k] = 0;
            continue ;
        }
        else
        {
            double sum = 0.0;
            double gama = 0.0;
            sum =                        x[i+1]* b[0][0][1] +   x[i-m+1]* b[0][0][2] +   x[i-m]* b[0][0][3]+   x[i-1]* b[0][0][4] +   x[i+m-1]* b[0][0][5] +   x[i+m]* b[0][0][6] +
                  + x[i+k]* b[0][1][0] + x[i+1+k]* b[0][1][1] + x[i-m+1+k]* b[0][1][2] + x[i-m+k]* b[0][1][3] + x[i-1+k]* b[0][1][4] + x[i+m-1+k]* b[0][1][5] + x[i+m+k]* b[0][1][6];
            gama = 0.0;
            if (column == 1)
            {
                if (row == 1)
                    gama = sum - f1*h/2;
                else if (row < n)
                    gama = sum - f1*h;
                else if (row == n)
                    gama = sum - f1*h/2.0;
            }
            else if (column == n)
            {
                if (row == 1)
                    gama = sum - f2*h/2;
                else if (row < n)
                    gama = sum - f2*h;
                else
                    gama = sum - f2*h/2;
            }
            else
            {
                gama = sum;
            }

            if (row == n && friction)
            {
                double gg = g_tr * h;
                if (column == 1 || column == n)
                {
                    gg = gg / 2;
                }
                if (gama <= gg)
                    gama += gg;
                else if (gama >= gg)
                    gama -= gg;
                else
                    gama = 0;
            }
            x_new = -(gama) / b[0][0][0];

            newdx = abs(x_new - x[i]);
            if (newdx > maxdx)
                maxdx = newdx;
            x[i] = x_new;

            //cout << "; x[i] = " << x[i];
            sum = x[i]*b[1][0][0] + x[i+1]*b[1][0][1] +   x[i-m+1]*b[1][0][2] +   x[i-m]*b[1][0][3] +   x[i-1]*b[1][0][4] +   x[i+m-1]*b[1][0][5] +   x[i+m]*b[1][0][6] +
                                  + x[i+1+k]*b[1][1][1] + x[i-m+1+k]*b[1][1][2] + x[i-m+k]*b[1][1][3] + x[i-1+k]*b[1][1][4] + x[i+m-1+k]*b[1][1][5] + x[i+m+k]*b[1][1][6];
            gama = sum;
            if (row == n && lagrange)
            {
               double temp = l[column] * h;
               if (column == 1 || column == n)
                   temp = temp / 2;
               gama = gama - temp;
            }
            y_new = -(gama) / b[1][1][0];

            if (row == n && (!lagrange) && unpenetration && y_new < 0)
                y_new = 0.0;


            newdx = abs(y_new - x[i+k]);
            if (newdx > maxdx)
                maxdx = newdx;
            x[i+k] = y_new;
        }

//        if (row == 1)

//        else
//        {
//            double sum = x[i]*b[1][0][0] + x[i+1]*b[1][0][1] +   x[i-m+1]*b[1][0][2] +   x[i-m]*b[1][0][3] +   x[i-1]*b[1][0][4] +   x[i+m-1]*b[1][0][5] +   x[i+m]*b[1][0][6] +
//                                         + x[i+1+k]*b[1][1][1] + x[i-m+1+k]*b[1][1][2] + x[i-m+k]*b[1][1][3] + x[i-1+k]*b[1][1][4] + x[i+m-1+k]*b[1][1][5] + x[i+m+k]*b[1][1][6];
//            double gamma = sum;
//            double x_new = -(gamma) / b[1][1][0];
//            if (unpenetration && x_new < 0)
//                x_new = 0;
//            double newdx = abs(x_new - x[i+k]);
//            if (newdx > maxdx)
//                maxdx = newdx;
//            x[i+k] = x_new;
//            //cout << "; < x[i+k] = " << x[i+k] << endl;
//        }
    }
    return maxdx;
}

void QuadGaussSeidelSolver::showGraph()
{
    QObject *par;
    par = new QObject();
    QProcess proc(par);

    proc.setWorkingDirectory(qApp->applicationDirPath() + QDir::separator());
    proc.start("python build-graph.py");
    proc.waitForStarted();
    proc.write(QString("%1\n").arg(2).toStdString().c_str());
    proc.write(QString("%1\n").arg(n).toStdString().c_str());
    for (int i = 0; i < 2*k; i++)
    {
        proc.waitForBytesWritten();
        proc.write(QString("%1\n").arg(x[i]).toStdString().c_str());
    }
    proc.waitForFinished(-1);

}

double QuadGaussSeidelSolver::updateLagrange(double ro)
{
    double outer_dx = 0.0;
    for (int i = 1; i < m-1; i++)
    {
        int ind = m*n+i;
        l[i] = max(0.0, l[i] - ro * x[ind+k]);
    }
    for (int i = 0; i < m; i++)
    {
        if (abs(l[i] - lOld[i]) > outer_dx)
            outer_dx = abs(l[i] - lOld[i]);
        lOld[i] = l[i];
    }
    return outer_dx;
}
